Magnet hypothesis: plant-pollinator interactions
Purpose: A test of the magnet hypothesis was examined in Mojave National Preserve by Ally Ruttan.
Hypothesis: Floral resource island created by shrubs and the associated beneficiary annual plants will positively and non-additively influence pollinator visitation rates.
Predictions:
(1) The frequency of pollinator visitations to annuals is greater under shrubs than in the paired-open microsites.
(2) Annual plants under flowering entomophilous shrubs (Larrea tridentata) will have a higher frequency of pollinator visitations than annual plants under anemophilous shrubs (Ambrosia dumosa) because of higher concentrations of floral resources for pollinators.
(3) Shrubs with annuals in their understory will have a higher frequency of pollinator visitations than shrubs without annuals due to increased concentrations of floral resources for pollinators.
(4) Sites with both shrubs and annuals will have the highest frequency of pollinator visitations to both the shrubs and the annuals.
An interesting corollary is that there are appropriate floral resources for desert pollinators, that they discriminate, and that entomophilous and anemophilous shrubs facilitate flowering similarly.
Data wrangling
#libraries
library(tidyverse)
library(DT)
library(lubridate)
#meta-data
meta <- read_csv("data/meta-data.csv")
datatable(meta)#data
data.2015 <- read_csv("data/MNP.2015.csv")
data.2016 <- read_csv("data/MNP.2016.csv")
#merge
data <- rbind.data.frame(data.2015, data.2016)
data <- data %>% rename(net.treatment = treatment) %>% na.omit(data)
data$year <- as.character(data$year)
data$rep <- as.character(data$rep)
#tidy data to expand treatment column (current structure is a mix of three factors)
#microsite
data <- data %>% mutate(microsite = ifelse(net.treatment %in% c("SA", "SAA", "SX"), "Larrea", ifelse(net.treatment %in% c("OA"), "open", ifelse(net.treatment %in% c("AMB"), "Ambrosia", NA))))
length(unique(data$microsite))## [1] 3#annuals
data <- data %>% mutate(annuals = ifelse(net.treatment %in% c("SA", "SAA"), "annuals", ifelse(net.treatment %in% c("OA"), "annuals", ifelse(net.treatment %in% c("AMB"), "annuals", "none"))))
length(unique(data$annuals))## [1] 2#flowering shrub
data <- data %>% mutate(flowering.shrub = ifelse(net.treatment %in% c("SAA", "SX"), "flowering shrub", ifelse(net.treatment %in% c("AMB", "SA", "OA"), "no shrub flowers", "NA")))
length(unique(data$annuals))## [1] 2#visitation duration needs to be weighted by duration (i.e. total recording time to address sampling effort)
data <- data %>% mutate(weighted.visitation.duration = as.duration(visitation.duration)/as.duration(total.duration))
#frequency counts in a separate dataframe (weighted by duration of recording)
frequency <- data %>% group_by(year, day, net.treatment, rep, microsite, annuals, flowering.shrub) %>% summarise(net.time = sum(total.duration), mean.visitation.duration = mean(weighted.visitation.duration), mean.floral.density = mean(floral.density), count = n())
frequency$net.time <- as.numeric(frequency$net.time) #converts to total seconds
frequency$rate <- as.numeric(frequency$count)/frequency$net.time #weight frequency by net time
datatable(frequency)Data visualization
#higher-order treatment patterns in frequency####
ggplot(frequency, aes(net.treatment, rate)) + geom_boxplot() + facet_wrap(~year)
ggplot(frequency, aes(microsite, rate)) + geom_boxplot() + facet_wrap(~year*annuals)
ggplot(frequency, aes(microsite, rate)) + geom_boxplot() + facet_wrap(~year*annuals*flowering.shrub)
#relationships with sampling effort
ggplot(frequency, aes(net.time, count, color = year)) + geom_point()
ggplot(frequency, aes(net.time, count, color = year)) + geom_point() + facet_wrap(~microsite)
#floral density
ggplot(frequency, aes(mean.floral.density, rate, color = year)) + geom_point()
ggplot(frequency, aes(mean.floral.density, rate, color = microsite)) + geom_point() + facet_wrap(~year)
ggplot(frequency, aes(mean.floral.density, rate, color = microsite)) + geom_point() + facet_wrap(~year)
#mean visitation duration is also really important because you have more visits or longer visits####
ggplot(frequency, aes(net.treatment, mean.visitation.duration)) + geom_boxplot() + facet_wrap(~year)
ggplot(frequency, aes(microsite, mean.visitation.duration)) + geom_boxplot() + facet_wrap(~year*annuals)
ggplot(frequency, aes(microsite, mean.visitation.duration)) + geom_boxplot() + facet_wrap(~year*annuals*flowering.shrub)
#relationships with sampling effort
ggplot(frequency, aes(net.time, mean.visitation.duration, color = year)) + geom_point()
ggplot(frequency, aes(net.time, mean.visitation.duration, color = year)) + geom_point() + facet_wrap(~microsite)
#floral density
ggplot(frequency, aes(mean.floral.density, mean.visitation.duration, color = year)) + geom_point()
ggplot(frequency, aes(mean.floral.density, mean.visitation.duration, color = microsite)) + geom_point() + facet_wrap(~year)
ggplot(frequency, aes(mean.floral.density, mean.visitation.duration, color = microsite)) + geom_point() + facet_wrap(~year)
EDA
#test distributions and explore outliers
summary(frequency)## year day net.treatment
## Length:266 Length:266 Length:266
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## rep microsite annuals
## Length:266 Length:266 Length:266
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## flowering.shrub net.time mean.visitation.duration
## Length:266 Min. : 300 Min. :0.0000000
## Class :character 1st Qu.: 4500 1st Qu.:0.0008142
## Mode :character Median : 23010 Median :0.0065703
## Mean : 64855 Mean :0.0159156
## 3rd Qu.:100660 3rd Qu.:0.0124122
## Max. :542929 Max. :0.4941065
## mean.floral.density count rate
## Min. : 5.00 Min. : 1.00 Min. :0.0001774
## 1st Qu.: 23.13 1st Qu.: 3.00 1st Qu.:0.0002130
## Median : 71.26 Median : 8.00 Median :0.0002780
## Mean : 98.74 Mean : 15.14 Mean :0.0005540
## 3rd Qu.:200.00 3rd Qu.: 22.00 3rd Qu.:0.0011111
## Max. :200.00 Max. :109.00 Max. :0.0033333#check orthogonality
freq.2015 <- frequency %>% filter(year == 2015)
freq.2016 <- frequency %>% filter(year == 2016)
#annual treatment
ggplot(freq.2015, aes(rate, fill = annuals)) + geom_histogram() + facet_wrap(~microsite)
#only larrea tested annuals and no-annuals in 2015
ggplot(freq.2016, aes(rate, fill = annuals)) + geom_histogram() + facet_wrap(~microsite)
#only larrea tested annuals and no-annuals in 2016
#flowering shrubs
ggplot(freq.2015, aes(rate, fill = flowering.shrub)) + geom_histogram() + facet_wrap(~microsite)
#ha, only larrea in 2015
ggplot(freq.2016, aes(rate, fill = flowering.shrub)) + geom_histogram() + facet_wrap(~microsite)
#ok, only larrea again but had ambrosia
#conclusion, separate years
require(fitdistrplus)
descdist(freq.2015$rate, boot = 1000)
## summary statistics
## ------
## min: 0.0001873882 max: 0.001337793
## median: 0.0003267998
## mean: 0.0004887686
## estimated sd: 0.0003435598
## estimated skewness: 1.110484
## estimated kurtosis: 2.810992descdist(freq.2016$rate, boot = 1000)
## summary statistics
## ------
## min: 0.0001773679 max: 0.003333333
## median: 0.0002480947
## mean: 0.0006144842
## estimated sd: 0.000504536
## estimated skewness: 1.331961
## estimated kurtosis: 7.146892detach("package:fitdistrplus", unload = TRUE)Models
#GLM for count and weight by net.time (alt - use MASS and glm.nb)
#library(MASS) #need for glm.nb
#all codes
#2015
m <- glm(count~net.treatment, family = "poisson", weight = net.time, data = freq.2015)
anova(m, test = "Chisq") ## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 127 47238655
## net.treatment 3 13736996 124 33501659 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ggplot(freq.2015, aes(count, fill = microsite, weight = net.time)) + geom_histogram(binwidth = 10) + scale_fill_brewer(palette = "Blues") #not a bad plot
#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")## $lsmeans
## net.treatment lsmean SE df asymp.LCL asymp.UCL
## OA 3.333027 0.0001166025 NA 3.332799 3.333256
## SA 2.118555 0.0005763230 NA 2.117425 2.119684
## SAA 2.993104 0.0001650586 NA 2.992780 2.993427
## SX 1.804488 0.0007044162 NA 1.803107 1.805869
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## OA - SA 1.2144727 0.0005880003 NA 2065.429 <.0001
## OA - SAA 0.3399236 0.0002020903 NA 1682.038 <.0001
## OA - SX 1.5285393 0.0007140016 NA 2140.806 <.0001
## SA - SAA -0.8745491 0.0005994936 NA -1458.813 <.0001
## SA - SX 0.3140666 0.0009101376 NA 345.076 <.0001
## SAA - SX 1.1886157 0.0007234961 NA 1642.878 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimates#2016
m <- glm(count~net.treatment, family = "poisson", weight = net.time, data = freq.2016)
anova(m, test = "Chisq") ## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 137 174945050
## net.treatment 4 15411680 133 159533369 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ggplot(freq.2016, aes(count, fill = microsite, weight = net.time)) + geom_histogram(binwidth = 10) + scale_fill_brewer(palette = "Blues") #not a bad plot
#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")## $lsmeans
## net.treatment lsmean SE df asymp.LCL asymp.UCL
## AMB 3.713592 8.577550e-05 NA 3.713424 3.713760
## OA 3.834869 7.552004e-05 NA 3.834721 3.835017
## SA 1.176321 2.278620e-03 NA 1.171855 1.180787
## SAA 3.970907 6.231829e-05 NA 3.970785 3.971029
## SX 1.432814 1.776673e-03 NA 1.429332 1.436297
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## AMB - OA -0.1212771 1.142835e-04 NA -1061.195 <.0001
## AMB - SA 2.5372705 2.280234e-03 NA 1112.724 <.0001
## AMB - SAA -0.2573154 1.060236e-04 NA -2426.963 <.0001
## AMB - SX 2.2807774 1.778742e-03 NA 1282.242 <.0001
## OA - SA 2.6585475 2.279871e-03 NA 1166.096 <.0001
## OA - SAA -0.1360383 9.791244e-05 NA -1389.388 <.0001
## OA - SX 2.4020544 1.778277e-03 NA 1350.776 <.0001
## SA - SAA -2.7945859 2.279472e-03 NA -1225.980 <.0001
## SA - SX -0.2564931 2.889407e-03 NA -88.770 <.0001
## SAA - SX 2.5380928 1.777765e-03 NA 1427.687 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates#treatments split out